Generalized-Ensemble Algorithms for the Isobaric-Isothermal Ensemble 
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We present generalized-ensemble algorithms for isobaric-isothermal molecular simulations. In 
addition to the multibaric-multithermal algorithm and replica-exchange method for the isobaric- 
isothermal ensemble, which have already been proposed, we propose a simulated tempering method 
for this ensemble. We performed molecular dynamics simulations with these algorithms for an 
alanine dipeptide system in explicit water molecules to test the effectiveness of the algorithms. 
We found that these generalized-ensemble algorithms are all useful for conformational sampling of 
biomolecular systems in the isobaric-isothermal ensemble. 

PACS numbers: 

Monte Carlo (MC) and molecular dynamics (MD) simulations with generalized-ensemble algorithms have been 
widely used for studies of proteins and peptides in order to have efficient conformational sampling (for a review, see, 
e.g., ref. 1). Three well-known generalized-ensemble algorithms are the multicanonical algorithm (MUCA) @, Q (the 
MD version was developed in refs. 4 and 5), replica-exchange method (REM) [|| (the method is also referred to as 
parallel tempering Ml and the MD version, which is referred to as REMD, was developed in ref. 8), and simulated 
tempering (ST) [9|, llCf . These methods can be used for obtaining accurate physical quantities in the canonical 
ensemble. Among them, REM (particularly REMD) is often used because the weight factor is a priori known (i.e., 
the Boltzmann factor), while those for MUCA and ST have to be determined before the simulations. 

Multidimensional (or multivariable) extensions of the original generalized-ensemble algorithms have been developed 
in many ways (see the references in ref. 1) and recently the general formulations have been given in refs. 11-13. 

In this letter, we present a multidimensional simulated tempering algorithm in the isobaric-isothermal (NPT) 
ensemble. Several generalized-ensemble algorithms for the NPT ensemble have been proposed (for example, the 
multibaric-multithermal (MUBATH) algorithm [l4]-[l7j and REM in temperature and pressure [IBt21|) and here, we 
present an ST algorithm for the NPT ensemble to complete generalized-ensemble algorithms for this ensemble. 

We first briefly review the generalized-ensemble algorithms for isobaric-isothermal molecular simulations. Let us 
consider a physical system that consists of N atoms and that is in a box of a finite volume V . The states of the 
system are specified by coordinates r = {r±, r%, ■ • ■ , r } and momenta p = {p\ , P2 , ■ • ■ , Pn } of the atoms and volume 
V of the box. The potential energy E(r, V) for the system is a function of r and V. 

We first describe MC simulation algorithms for MUCA, REM, and ST in the NPT ensemble. In these cases, 
momenta of atoms do not have to be considered. To make a system an equilibrium state, the detailed balance 
condition is imposed and a transition probability w(X —> X') from an old state X to a new state X' can be given by 
the Metropolis criterion. [22j 

In MUBATH simulations, we introduce a function fi(E, V) and use a weight factor W m bt(E, V) = exp [—/3oH(E, V)] 
so that the distribution function f m bt{E, V) of E and V may be uniform: 

f mht (E, V) oc n(E, V)W mht {E, V) = constant, (1) 

where /3q is an arbitrary inverse reference temperature defined as /3q — I/ZcbTo (fee is the Boltzmann constant) and 
n(E, V) is the density of states. 

To perform MUBATH MC simulations, the trial moves are generated in the same way as in the usual constant 
NPT MC simulations HU and the transition probability from X = {s, V} to X' = {s', V'} is given by QUI 

w mht (X -> X') = min [l,exp(-Anbt)] , (2) 

where 

A„bt = Po {U [E(s', V), V] - U [E(s, V), V] - Nk B T HV'/V)} , (3) 

and s = {si, S2, • ■ • , s^} is the scaled coordinates defined by Sj = V -1 / 3 ^ [i = 1, 2, • • • , N). Here, we are assuming 
the box is a cube of side U -1 / 3 . 

In REM simulations, we prepare a system that consists of Mt x Mp non-interacting replicas of the original system, 
where Mt and Mp are the number of temperature and pressure values used in the simulation, respectively. The 
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replicas are specified by labels i (i — 1, 2, ■ ■ • , Mt x Mp), temperature by m t (m t = 1, 2, ■ ■ • , Mt), and pressure by 
m p (m p = 1,2, ••• ,M P ). 

To perform REM MC (REMC) simulations, we carry out the following two steps alternately: (1) perform a 
usual constant NPT simulation in each replica at assigned temperature and pressure and (2) try to exchange the 
replicas. If the temperature (specified by m t and n t ) and pressure (specified by m p and n p ) between the repli- 
cas are exchanged, the transition probability from X = {••• , (s^, VM; T mt , P OT ), • ■ • ,(s^\V^;T nt ,P n ),■••} to 
X' = {■■-, (aM , FW ; T nt , P„„), • • • , (sM , ; T mt , P m „), • • • } at the trial is given by [H Q 

w rom (A -)• X') = min [1, exp(-A-om)] , (4) 

where 

Aem = (Am - Pn t ) [E( S ^,V^) - E(s®, V W)] + (/3„ lt P„ lp - p nt P np ) (yW - V®) . (5) 

In ST simulations, we introduce a function g(T, P) and use a weight factor W s t(E, V; T, P) = exp[— j3(E + PV) + 
g(T, P)] so that the distribution function f s t(T, P) of T and P may be uniform: 

f st (T,P) cx f dV f dr W st [E(r,V),V;T,P] = constant. (6) 
Jo Jv 

From eq. (J6j , it is found that g(T,P) is formally given by 

g(T,P) = -1a^J dV J drexp[-(3{E{r,V) + PV)]j, (?) 

and the function is the dimensionless Gibbs free energy except for a constant. 

To perform ST MC simulations, we carry out the following two steps alternately: (1) perform a usual constant NPT 
simulation and (2) try to update the temperature and pressure. The transition probability from X = {s, V; T, P} to 
X 1 = {s, V; T', P'} for this trial is given by 

tu rt (A- mm [l,exp(- Art)], (8) 

where 

At = W - P)E(s, V) + (P'P 1 - PP)V - [g(r, P>) - g(T, P)] . (9) 

For MD simulations with MUCA, REM, and ST in the NPT ensemble, the actual formulations depend on constant 
temperature and pressure algorithms. Here, we employ the MD methods with the Martyna- Tobias-Klein (MTK) 
algorithm, [24j whose equations of motion follow Nose [1^, [26| and Hoover (27} for the thermostat and Andersen (28[ 
for the barostat. 

To perform MUBATH MD simulations, we solve the usual equations of motion for the MTK algorithm except that 
it is necessary to modify the equations for pi and p £ as follows: 

d Pl dUdE ( 1\ Pe p t 



n 



" 3y -dv + N^^-Q p - (10b) 
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where 

P 1 



"tf_dH(" r _dE d_E' 
^ nn dE \ 1 dn dV 

,i=l \i=l / 



(11) 



p £ is the momentum associated with the logarithm of V, p^ is the momentum of the thermostat, and W and Q are 
the masses of barostat and thermostat, respectively. 

When we perform MD simulations with REM and ST, the momenta should be rescaled if the replicas are exchanged 
for the temperature in REM Q and the temperature is updated in ST [12j . 

In REM MD (REMD) simulations with the MTK algorithm, let an old state be A = 
{•■• ,{rW,V®,pW,p®,pf;T mt ,P mp ),.-- ,(r^,V^,p^,p^ ] ,pf;T nt ,P np ),---} and a new state X = 
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{•■• ,(r®,VW,pW,p®',pf';T nt ,P np ),--- ,(r^,V^,p^',p^',pf;T mt ,P mp ),---}. When replica exchange is 
performed, the momenta should be rescaled as follows: [291 ] 



T 




/ T nt W n r-l r-, T mt W m r-i 

4 = V^: 4 ' 4 = v^r pi (12b) 



p M/ = / ^Qn p f f W = / ^"LpW , (12c) 

where W m and <5 m are the masses in the simulation at T mt and P mp , and W„ and Q n are the ones at T nt and P„ p . 
Then the transition probability (Q} and (J5J) can also be used in the REMD simulations. 

In ST MD simulations with the MTK algorithm, let an old state be X = {r,V,p,p e ,p^;T, P} and a new state 
X' = {r, V,p',p' s ,p'^;T', P'}. lfp',p' e and p'^ are written as 



, T> , T'W , T'Q' 

where W and Q' are the masses in the simulation at T 1 and P', then it can be shown that the transition probability 
is given by eqs. ([8]) and (|9]). 

While the weight factors in REM simulations are a priori known, those in MUBATH and ST simulations have to 
be determined before the simulations. We propose two-dimensional extensions of the replica-exchange multicanonical 
algorithm (REMUCA) [3(| and the replica-exchange simulated tempering (REST) [3l| for the isobaric-isothermal 
ensemble to overcome the difficulty of determining the weight factors. First, we perform a REM simulation for the 
isobaric-isothermal ensemble and then calculate the density of states and the dimensionless Gibbs free energy by the 
multiple-histogram reweighting techniques, [Hj] or the weighted histogram analysis method (WHAM), [33| using the 
results of the REM simulation. The WHAM equations for isobaric-isothermal simulations are given by 

M T M P 

Y / Y j (l+2r ij )- 1 N ij (E,V) 

n ^ V ) - TI^^ ■ ( 14a ) 

2 J2 + exp [ 9ij - 0t(E + PjV)] 

i=l j=l 

exp {- 9ij ) =Y,Y, V ^ ex P {-&( E + P ^ ' ( 14b ) 

E V 

where gij = g(Ti, Pj), Nij(E, V) is the histogram of E and V, Uij is the number of the data, and Tjj is the integrated 
autocorrelation time in the simulation at Tj and Pj , respectively. These equations need to be solved self-consistcntly. 
We remark that for biomolecular systems, it was found that may be set to a constant value for all i and j in eq. 

(Ha) . 13 

In the isobaric-isothermal version of REST, the weight factor can also be determined using the dimensionless Gibbs 
free energy obtained from eqs. (|14a[) and (|14b[) . 

In order to verify that the generalized-ensemble algorithms discussed above can be effective for conformational 
sampling and give the same results, we performed MD simulations with the three generalized-ensemble algorithms. 
We used a system of an alanine dipeptide in 73 surrounding water molecules and the system was placed in a cubic 
cell with periodic boundary conditions. Both of the backbone dihedral angles <f> and ip of the peptide were initially 
set to 180°. The simulations were performed with the tinker program package. [34[ Several of the programs were 
modified and a few programs were added so that MUBATH, REM, and ST simulations with the MTK algorithm can 
be performed. We used the AMBER parm99 force field [3|| for the peptide and the TIP3P water model. [36[ The 
electrostatic interactions were calculated by the particle mesh Ewald method. [37], [38j] In the van der Waals interaction 
calculations, we used the spherical cutoff method and the cutoff distance was set to 12.0 A. The integration of the 
equations of motion was employed by the method proposed by Martyna et al. [39[ and the Shake/Rattle/Roll constraint 
method [39M4l| was used so that the water molecules are rigid body molecules. The unit time step was set to 0.5 fs. 
The mass parameters W and Q were determined as in ref. 39. 
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FIG. 1: Results of the two-dimensional REMD simulation: (a) the time series of T and P in one of the replicas and (b) the 
time series of the replica label at 280 K and 250 MPa. 



First, we performed the two-dimensional REMD simulation. The simulation time was set to 2.0 ns. We used 
the following 6 temperature (Ti, • ■ • ,T 6 ) and 4 pressure (Pi, • • • , P 4 ) values: 280, 305, 332, 362, 395, and 430 K for 
temperature and 0.1, 65, 150, and 250 MPa for pressure. At the replica-exchange trial, either exchanging temperature 
(T-exchange) or exchanging pressure (P-exchange) was chosen randomly and then the pairs {(Ti, T2), (X3, T4), (T5, Tq)} 
or {(T 2 , T3), (T4, T5)} for T-exchange and {(Pi, P2), (P3, Pa)} or {(P2, P3)} for P-exchange were also chosen randomly. 

Figure [lja) shows the time series of T and P in one of the replicas and Fig. [T{b) shows the time series of the label 
of the replicas in the simulation at 280 K and 250 MPa. From these figures, it is found that random walks in T-P 
space were realized in each of the replicas. The acceptance ratios in the REM simulation were in the range from 
0.152 to 0.205 for T-exchange and from 0.296 to 0.446 for P-exchange, respectively, which also implies that sufficient 
number of replica exchange was realized. 

We then performed 24 MUBATH simulations of 2.0 ns, where the total simulation time was 48 ns so that it is equal 
to that in the REMD simulation. In each of the simulations, different initial velocities were given. The data were 
sampled every 100 fs. The reference temperature was set to 430 K. dTL/dE and dTL/dV were interpolated by the 
third-order polynomial, following ref. 17. 

Figure [5] shows the probability distributions of E and V from the MUBATH simulations. From Fig. [Ha), it is found 
that the MUBATH simulations gave a uniform distribution in the range where the density of states was obtained 
accurately in the REMD simulation and that the distribution of the MUBATH simulations was much wider than the 
ones for the NPT ensemble (see Fig. [2fb)). 

Twenty-four ST simulations of 2.0 ns were carried out, which gave us the same number of sampled data as in 
the REMD and MUBATH simulations. In each simulation, different initial velocities were given. We used the same 
temperature and pressure values as in the REMD simulation. We tried to update the parameter (T or P) every 100 
fs and the data were sampled just before the trials. The choice between updating T (T-update) and updating P 
(P-update) was made randomly and then either T mt -i or T mt+ i for T-update and P m _i or P m + i for P-update 
was also chosen randomly, where mt and m p stand for the labels corresponding to the temperature and the pressure 
before the trial, respectively. 




FIG. 2: Logarithm of probability distributions / of E and V (a) in the MUBATH simulations and (b) for the NPT ensemble 
at 280 K and 250 MPa (left) and at 430 K and 0.1 MPa (right). The distributions corresponding to the NPT ensemble were 
calculated by the single-histogram reweighting techniques [42| using the results of the MUBATH simulations. 
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FIG. 3: Results of the ST simulations: (a) the time series of T and P for 2.0 ns and (b) the logarithm of the probability 
distribution / of T and P. 



Figure 03a) shows the time series of T and P and Fig. [3Jb) shows the probability distribution of T and P. These 
figures indicate that random walks in T-P space were successfully realized. The acceptance ratios in the ST simulations 
were in the range from 0.260 to 0.430 for T-update and from 0.426 to 0.642 for P-update, respectively. 

Figure 0] shows the probability distributions of the backbone dihedral angles at 298 K and 0.1 MPa in all the 
simulations. Compared with the simulations with REMD, MUBATH, and ST, we also performed 24 conventional 
isobaric-isothermal simulations of 2.0 ns at 298 K and 0.1 MPa with different initial velocities. The dihedral angle 
distributions in the simulations with the generalized-ensemble algorithms had a small peak in 0° < <j> < 90° and 
—90° < ip < 90°, although there was no peak in the range in the distribution of the conventional simulations. All the 
simulations with the three generalized-ensemble algorithms were able to provide the same results and reproduce the 
distribution obtained previously in refs. 43 and 44. 

In this letter, we presented the extensions of MUCA, REM, and ST for isobaric-isothermal molecular simulations. 
These algorithms can be effective methods for conformational sampling and give accurate physical quantities in the 
isobaric-isothermal ensemble. Therefore, one can use the generalized-ensemble algorithms to study temperature and 
pressure effects on complex systems such as biomolecular systems. 

We expect the simulated tempering to be a more useful method for simulations in large systems. Many replicas are 
needed in REM simulations of large systems, while only one replica is used in ST simulations. In addition, ST can be 
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FIG. 4: (Color online) Contour maps of probability distribution of the backbone dihedral angles <j> and ip in the simulations 
with (a) REMD, (b) MUBATH, and (c) ST and (d) in the conventional isobaric-isothermal simulations. In these figures, the 
probability distributions at 298 K and 0.1 MPa are plotted in logarithmic scale and were calculated by the single-histogram 
reweighting techniques for MUBATH and by WHAM for REMD and ST. 
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easily implemented into existing program packages compared to MUCA (or MUBATH), in which a force calculation 
part of programs has to be modified for MD simulations. Therefore, molecular simulations with ST will be more 
easily applicable to larger systems than with the other two generalized-ensemble algorithms. 
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